library(rstanarm)
Loading required package: Rcpp
This is rstanarm version 2.21.1
- See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
- Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
- For execution on a local, multicore CPU with excess RAM we recommend calling
  options(mc.cores = parallel::detectCores())

Attaching package: ‘rstanarm’

The following object is masked from ‘package:rstan’:

    loo

lq_theta_y = function(sigmaSQ, beta0, beta1, lpr = df$lprice, lcar = df$lcarat) {
  ld = dnorm(beta0, 0, 100, log = TRUE) + 
    dnorm(beta1, 0, 100, log = TRUE) + 
    dinvgamma(sigmaSQ, shape = 0.01, rate = 0.01, log = TRUE) +
    sum(dnorm(lpr, mean = (beta0 + beta1 * lcar), sd = sqrt(sigmaSQ), log = TRUE))
  return(ld)
}

mh = function(B, theta_start) {
  theta = array(0, c((B+1), 3), dimnames = list(c(), c("sigmaSQ", "beta0", "beta1")))
  
  theta[1,1] = theta_start[1]
  theta[1,2] = theta_start[2]
  theta[1,3] = theta_start[3]

  for (i in 2:dim(theta)[1]) {
    
    ###  step for sigmaSQ
    beta0_star = theta[(i-1),2]
    beta1_star = theta[(i-1),3]
    
    sigmaSQ_star = rtruncnorm(n = 1, a=0, b=Inf, mean = theta[(i-1),1], sd = 0.1)

    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      log(dtruncnorm(x = sigmaSQ_star, a=0, b=Inf, mean = theta[(i-1),1], sd = 0.1))
    den_logr = lq_theta_y(sigmaSQ = theta[i-1, 1], beta0 = beta0_star, beta1 = beta1_star) -
      log(dtruncnorm(x = theta[i-1, 1], a=0, b=Inf, mean = sigmaSQ_star, sd = 0.1)) 
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,1] = sigmaSQ_star
    } else {
      theta[i,1] = theta[(i - 1), 1]
    }
    
    ###  step for beta0
    beta1_star = theta[(i-1),3] # it is the repeated code, but I need it to keep the interpretability
    sigmaSQ_star = theta[i,1] # update sigmaSQ_star after the Gibbs step for sigmaSQ
    
    beta0_star = rnorm(1, theta[(i-1),2], 0.1) # !!! check the parametrization (0.1 or 0.1^2)
    
    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      dnorm(x = beta0_star, mean = theta[(i-1),2], sd = 0.1, log = TRUE)
    den_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = theta[i-1, 2], beta1 = beta1_star) -
      dnorm(x = theta[(i-1),2], mean = beta0_star, sd = 0.1, log = TRUE)
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,2] = beta0_star
    } else {
      theta[i,2] = theta[(i - 1), 2]
    }
    
    ###  step for beta1
    sigmaSQ_star = theta[i,1] # it is the repeated code, but I need it to keep the interpretability
    beta0_star = theta[i,2] # update beta0_star after the Gibbs step for beta0

    beta1_star = rnorm(1, theta[(i-1),3], 0.1) # !!! check the parametrization (0.1 or 0.1^2)
    
    num_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = beta1_star) -
      dnorm(x = beta1_star, mean = theta[(i-1),3], sd = 0.1, log = TRUE)
    den_logr = lq_theta_y(sigmaSQ = sigmaSQ_star, beta0 = beta0_star, beta1 = theta[i-1, 3]) -
      dnorm(x = theta[(i-1),3], mean = beta1_star, sd = 0.1, log = TRUE)
## !!! is the mean correct? Should  we use in den_logr mean = sigmaSQ_star or mean = theta[i-1, 1]? 
# check c-componentwise-mh v1 line 68, 72
    logr = num_logr - den_logr
    if (log(runif(1)) <= min(logr, 0)) {
      theta[i,3] = beta1_star
    } else {
      theta[i,3] = theta[(i - 1), 3]
    }
    
  }
  return(theta)
}

B = 10^5 
keep = (B/2 + 1):(B + 1)
chain1 = mh(B, theta_start = c(0.1, -1, -1))
chain2 = mh(B, theta_start = c(0.3, 0, 0))
chain3 = mh(B, theta_start = c(0.5, -1, 1))
chain4 = mh(B, theta_start = c(0.2, 1, 1))

mc = mcmc.list(mcmc(chain1[keep,]), mcmc(chain2[keep,]),
                 mcmc(chain3[keep,]), mcmc(chain4[keep,]))
summary(mc)

Iterations = 1:50001
Thinning interval = 1 
Number of chains = 4 
Sample size per chain = 50001 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean      SD  Naive SE Time-series SE
sigmaSQ 0.1152 0.00925 2.068e-05      6.977e-05
beta0   8.3798 0.02240 5.008e-05      1.647e-04
beta1   1.5064 0.03135 7.011e-05      2.141e-04

2. Quantiles for each variable:

           2.5%    25%    50%    75%  97.5%
sigmaSQ 0.09853 0.1088 0.1145 0.1211 0.1348
beta0   8.33546 8.3649 8.3797 8.3948 8.4234
beta1   1.44522 1.4854 1.5064 1.5276 1.5680
gelman.diag(mc, autoburnin = FALSE)
Potential scale reduction factors:

        Point est. Upper C.I.
sigmaSQ       1.03       1.08
beta0         1.00       1.01
beta1         1.01       1.03

Multivariate psrf

1.03
# extracting the MCMC samples from the soda_fit object
samples = extract(r_fit_A)
ncycles = length(samples[[1]])
T <- length(df$price)

# each row of yrep is a sample from the pp distribution
yrep = matrix(0, ncol = T, nrow = ncycles)
for (i in seq_len(T)) {
  mui = samples$beta0 + samples$beta1 * df$lcarat[i]
  yrep[, i] = rnorm(ncycles, mean = mui, sd = sqrt(samples$sigmaSQ))
}

# approximate posterior predictive density for an observation
# with the same covariates as observation 1
#plot(density(yrep[,1]))
color_scheme_set("red")

# posterior predictive check
y = df$lprice
ppc_hist(y, yrep[sample(1:ncycles, 8),]) + ggtitle("Model A: histogram comparing the response to 8 replicated data sets")



# scatterplot of y vs average yrep
ppc_scatter_avg(y, yrep) + ggtitle("Model A: scatterplot comparing the response to the average replicated data set")


ppc_dens_overlay(y, yrep = yrep[sample(1:ncycles, 100),]) + ggtitle("Model A: density plot comparing the density of the response to the densities of 100 replicated data sets")


ppc_error_scatter_avg(y, yrep = yrep) + ggtitle("Model A: scatter plot comparing the response to the average predictive error")


ppc_intervals(y, yrep = yrep)


ppc_error_scatter_avg_vs_x(y = y, yrep = yrep, x = df$lcarat)

ppc_intervals(y = y, yrep = yrep, x = df$lcarat) + ggtitle("Model A: Posterior predictive intervals to the observed data values") + xlab("log carat") + ylab("log price")

# extracting the MCMC samples from the soda_fit object
samples = extract(r_fit_B)
ncycles = length(samples[[1]])
T <- length(df$price)

# each row of yrep is a sample from the pp distribution
yrep = matrix(0, ncol = T, nrow = ncycles)
for (i in seq_len(T)) {
  mui = samples$beta0 + samples$beta1 * df$carat[i] + samples$beta2 * df$carat[i] * df$carat[i]
  yrep[, i] = rnorm(ncycles, mean = mui, sd = sqrt(samples$sigmaSQ))
}

# approximate posterior predictive density for an observation
# with the same covariates as observation 1
#plot(density(yrep[,1]))
color_scheme_set("green")

# posterior predictive check
y = df$price
ppc_hist(y, yrep[sample(1:ncycles, 8),]) + ggtitle("Model B: histogram comparing the response to 8 replicated data sets")



# scatterplot of y vs average yrep
ppc_scatter_avg(y, yrep) + ggtitle("Model B: scatterplot comparing the response to the average replicated data set")


ppc_dens_overlay(y, yrep = yrep[sample(1:ncycles, 100),]) + ggtitle("Model B: density plot comparing the density of the response to the densities of 100 replicated data sets")


ppc_error_scatter_avg(y, yrep = yrep) + ggtitle("Model B: scatter plot comparing the response to the average predictive error")


ppc_intervals(y = y, yrep = yrep, x = df$carat) + ggtitle("Model B: Posterior predictive intervals to the observed data values") + xlab("carat") + ylab("price")

samples = extract(r_fit_D)
ncycles = length(samples[[1]])
T <- length(df$price)
df$notic = as.integer(df$inclusions == "noticeable")

# each row of yrep is a sample from the pp distribution
yrep = matrix(0, ncol = T, nrow = ncycles)
for (i in seq_len(T)) {
  mui = samples$beta0 + samples$beta1 * df$lcarat[i] + samples$beta2 * df$notic[i]
  yrep[, i] = rlnorm(ncycles, mean = mui, sd = sqrt(samples$sigmaSQ))
}

color_scheme_set("blue")
ppc_intervals_grouped(y = df$price, yrep = yrep, x = df$carat, prob = 0.5, prob_outer = .95, group = df$notic) + ggtitle("Model D: Posterior predictive intervals to the observed data values") + xlab("lcarat") + ylab("price")


find_ep = function(yrep, min = 0.1, max = 2.4, step=0.1) {
  ep = array(0, dim=c((max-min+step)/step, 5), dimnames = list(c(), c("carat_min", "carat_max", "2.5%","Mean","97.5%")))
  n=0
  for (i in seq(min+step, max, by=step)) {
    n=n+1
    subset = ((df$carat>(i-step)) & (df$carat<=i))
    yrep_subset = yrep[,subset]
    ep[n,1] = i-step
    ep[n,2] = i
    q = quantile(yrep_subset, probs = c(0.025, 0.975))
    ep[n,3] = q[[1]]
    ep[n,4] = mean(yrep_subset)
    ep[n,5] = q[[2]]
  }
return(ep)
}

find_ep(yrep)
temp = ppc_intervals_data(y = df$price, yrep = yrep, x = df$lcarat, prob = 0.5, prob_outer = .95, group = df$notic)
cat("STAN Model E WAIC", waic_loo(r_fit_E)[1])
cat("STAN Model E LOOIC", waic_loo(r_fit_E)[2])
samples = extract(r_fit_E)
ncycles = length(samples[[1]])
T <- length(df_covid$bs)

# each row of yrep is a sample from the pp distribution
yrep = matrix(0, ncol = T, nrow = ncycles)
for (i in seq_len(T)) {
  log_l = log(df_covid$population[i]) +
    samples$beta0 + 
    samples$beta1 * df_covid$income[i] + 
    samples$beta2 * df_covid$bs[i]
  yrep[, i] = rpois(ncycles, lambda = exp(log_l))
}

ppc_intervals(y = y, yrep = yrep, x = df_covid$state_abb) + ggtitle("Model A: Posterior predictive intervals to the observed data values") + xlab("State") + ylab("Covid deaths")
LS0tCnRpdGxlOiAiNzM5MyBFeGFtIDIgUiBDb2RlIEFuZHJlaSBNYXR2ZWV2IgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCmBgYHtyIERhdGEgYW5kIHNldCB1cH0KCmxpYnJhcnkoZGF0YS50YWJsZSkKbGlicmFyeSh0cnVuY25vcm0pCmxpYnJhcnkoY29kYSkKZGV2dG9vbHM6Omluc3RhbGxfZ2l0aHViKCJzdGFuLWRldi9sb28iKQpsaWJyYXJ5KGxvbykKbGlicmFyeShyc3RhbikgCmxpYnJhcnkoc2hpbnlzdGFuKQpsaWJyYXJ5KGludmdhbW1hKQpsaWJyYXJ5KGJheWVzcGxvdCkKbGlicmFyeShwb3N0ZXJpb3IpCmxpYnJhcnkocnN0YW5hcm0pCgoKZGV2dG9vbHM6Omluc3RhbGxfZ2l0aHViKCJqZnJlbmNoL2JheWVzdXRpbHMiKQoKc2V0d2QoJy9Vc2Vycy9BTS9Eb2N1bWVudHMvX0NVIE1hc3RlcnMvMjAyMCBmYWxsIEJheWVzaWFuXzczOTMvY29kZS9CYXllc2lhbl9TdGF0aXN0aWNzX0NsYXNzX0NvZGUvRXhhbV9jb2RlJykKI2RmMSA9IGxvYWQoImRpYW1vbmRzX3NpbXBsZS5yZGEiKQojcm0oZGYxLCBkaWFtb25kc19zaW1wbGUpCmRmID0gYmF5ZXN1dGlsczo6ZGlhbW9uZHNfc2ltcGxlCgpkZiRscHJpY2UgPSBsb2coZGYkcHJpY2UpCmRmJGxjYXJhdCA9IGxvZyhkZiRjYXJhdCkKCmRmX2NvdmlkID0gYmF5ZXN1dGlsczo6Y292aWRfZGVjNAoKYGBgCgpgYGB7ciBQcm9ibGVtIDEgbG9nIGRlbnNpdHkgRnVuY3Rpb259CgpscV90aGV0YV95ID0gZnVuY3Rpb24oc2lnbWFTUSwgYmV0YTAsIGJldGExLCBscHIgPSBkZiRscHJpY2UsIGxjYXIgPSBkZiRsY2FyYXQpIHsKICBsZCA9IGRub3JtKGJldGEwLCAwLCAxMDAsIGxvZyA9IFRSVUUpICsgCiAgICBkbm9ybShiZXRhMSwgMCwgMTAwLCBsb2cgPSBUUlVFKSArIAogICAgZGludmdhbW1hKHNpZ21hU1EsIHNoYXBlID0gMC4wMSwgcmF0ZSA9IDAuMDEsIGxvZyA9IFRSVUUpICsKICAgIHN1bShkbm9ybShscHIsIG1lYW4gPSAoYmV0YTAgKyBiZXRhMSAqIGxjYXIpLCBzZCA9IHNxcnQoc2lnbWFTUSksIGxvZyA9IFRSVUUpKQogIHJldHVybihsZCkKfQoKYGBgCgpgYGB7ciBQcm9ibGVtIDEgbWggRnVuY3Rpb259CgptaCA9IGZ1bmN0aW9uKEIsIHRoZXRhX3N0YXJ0KSB7CiAgdGhldGEgPSBhcnJheSgwLCBjKChCKzEpLCAzKSwgZGltbmFtZXMgPSBsaXN0KGMoKSwgYygic2lnbWFTUSIsICJiZXRhMCIsICJiZXRhMSIpKSkKICAKICB0aGV0YVsxLDFdID0gdGhldGFfc3RhcnRbMV0KICB0aGV0YVsxLDJdID0gdGhldGFfc3RhcnRbMl0KICB0aGV0YVsxLDNdID0gdGhldGFfc3RhcnRbM10KCiAgZm9yIChpIGluIDI6ZGltKHRoZXRhKVsxXSkgewogICAgCiAgICAjIyMgIHN0ZXAgZm9yIHNpZ21hU1EKICAgIGJldGEwX3N0YXIgPSB0aGV0YVsoaS0xKSwyXQogICAgYmV0YTFfc3RhciA9IHRoZXRhWyhpLTEpLDNdCiAgICAKICAgIHNpZ21hU1Ffc3RhciA9IHJ0cnVuY25vcm0obiA9IDEsIGE9MCwgYj1JbmYsIG1lYW4gPSB0aGV0YVsoaS0xKSwxXSwgc2QgPSAwLjEpCgogICAgbnVtX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSBzaWdtYVNRX3N0YXIsIGJldGEwID0gYmV0YTBfc3RhciwgYmV0YTEgPSBiZXRhMV9zdGFyKSAtCiAgICAgIGxvZyhkdHJ1bmNub3JtKHggPSBzaWdtYVNRX3N0YXIsIGE9MCwgYj1JbmYsIG1lYW4gPSB0aGV0YVsoaS0xKSwxXSwgc2QgPSAwLjEpKQogICAgZGVuX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSB0aGV0YVtpLTEsIDFdLCBiZXRhMCA9IGJldGEwX3N0YXIsIGJldGExID0gYmV0YTFfc3RhcikgLQogICAgICBsb2coZHRydW5jbm9ybSh4ID0gdGhldGFbaS0xLCAxXSwgYT0wLCBiPUluZiwgbWVhbiA9IHNpZ21hU1Ffc3Rhciwgc2QgPSAwLjEpKSAKIyMgISEhIGlzIHRoZSBtZWFuIGNvcnJlY3Q/IFNob3VsZCAgd2UgdXNlIGluIGRlbl9sb2dyIG1lYW4gPSBzaWdtYVNRX3N0YXIgb3IgbWVhbiA9IHRoZXRhW2ktMSwgMV0/IAojIGNoZWNrIGMtY29tcG9uZW50d2lzZS1taCB2MSBsaW5lIDY4LCA3MgogICAgbG9nciA9IG51bV9sb2dyIC0gZGVuX2xvZ3IKICAgIGlmIChsb2cocnVuaWYoMSkpIDw9IG1pbihsb2dyLCAwKSkgewogICAgICB0aGV0YVtpLDFdID0gc2lnbWFTUV9zdGFyCiAgICB9IGVsc2UgewogICAgICB0aGV0YVtpLDFdID0gdGhldGFbKGkgLSAxKSwgMV0KICAgIH0KICAgIAogICAgIyMjICBzdGVwIGZvciBiZXRhMAogICAgYmV0YTFfc3RhciA9IHRoZXRhWyhpLTEpLDNdICMgaXQgaXMgdGhlIHJlcGVhdGVkIGNvZGUsIGJ1dCBJIG5lZWQgaXQgdG8ga2VlcCB0aGUgaW50ZXJwcmV0YWJpbGl0eQogICAgc2lnbWFTUV9zdGFyID0gdGhldGFbaSwxXSAjIHVwZGF0ZSBzaWdtYVNRX3N0YXIgYWZ0ZXIgdGhlIEdpYmJzIHN0ZXAgZm9yIHNpZ21hU1EKICAgIAogICAgYmV0YTBfc3RhciA9IHJub3JtKDEsIHRoZXRhWyhpLTEpLDJdLCAwLjEpICMgISEhIGNoZWNrIHRoZSBwYXJhbWV0cml6YXRpb24gKDAuMSBvciAwLjFeMikKICAgIAogICAgbnVtX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSBzaWdtYVNRX3N0YXIsIGJldGEwID0gYmV0YTBfc3RhciwgYmV0YTEgPSBiZXRhMV9zdGFyKSAtCiAgICAgIGRub3JtKHggPSBiZXRhMF9zdGFyLCBtZWFuID0gdGhldGFbKGktMSksMl0sIHNkID0gMC4xLCBsb2cgPSBUUlVFKQogICAgZGVuX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSBzaWdtYVNRX3N0YXIsIGJldGEwID0gdGhldGFbaS0xLCAyXSwgYmV0YTEgPSBiZXRhMV9zdGFyKSAtCiAgICAgIGRub3JtKHggPSB0aGV0YVsoaS0xKSwyXSwgbWVhbiA9IGJldGEwX3N0YXIsIHNkID0gMC4xLCBsb2cgPSBUUlVFKQojIyAhISEgaXMgdGhlIG1lYW4gY29ycmVjdD8gU2hvdWxkICB3ZSB1c2UgaW4gZGVuX2xvZ3IgbWVhbiA9IHNpZ21hU1Ffc3RhciBvciBtZWFuID0gdGhldGFbaS0xLCAxXT8gCiMgY2hlY2sgYy1jb21wb25lbnR3aXNlLW1oIHYxIGxpbmUgNjgsIDcyCiAgICBsb2dyID0gbnVtX2xvZ3IgLSBkZW5fbG9ncgogICAgaWYgKGxvZyhydW5pZigxKSkgPD0gbWluKGxvZ3IsIDApKSB7CiAgICAgIHRoZXRhW2ksMl0gPSBiZXRhMF9zdGFyCiAgICB9IGVsc2UgewogICAgICB0aGV0YVtpLDJdID0gdGhldGFbKGkgLSAxKSwgMl0KICAgIH0KICAgIAogICAgIyMjICBzdGVwIGZvciBiZXRhMQogICAgc2lnbWFTUV9zdGFyID0gdGhldGFbaSwxXSAjIGl0IGlzIHRoZSByZXBlYXRlZCBjb2RlLCBidXQgSSBuZWVkIGl0IHRvIGtlZXAgdGhlIGludGVycHJldGFiaWxpdHkKICAgIGJldGEwX3N0YXIgPSB0aGV0YVtpLDJdICMgdXBkYXRlIGJldGEwX3N0YXIgYWZ0ZXIgdGhlIEdpYmJzIHN0ZXAgZm9yIGJldGEwCgogICAgYmV0YTFfc3RhciA9IHJub3JtKDEsIHRoZXRhWyhpLTEpLDNdLCAwLjEpICMgISEhIGNoZWNrIHRoZSBwYXJhbWV0cml6YXRpb24gKDAuMSBvciAwLjFeMikKICAgIAogICAgbnVtX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSBzaWdtYVNRX3N0YXIsIGJldGEwID0gYmV0YTBfc3RhciwgYmV0YTEgPSBiZXRhMV9zdGFyKSAtCiAgICAgIGRub3JtKHggPSBiZXRhMV9zdGFyLCBtZWFuID0gdGhldGFbKGktMSksM10sIHNkID0gMC4xLCBsb2cgPSBUUlVFKQogICAgZGVuX2xvZ3IgPSBscV90aGV0YV95KHNpZ21hU1EgPSBzaWdtYVNRX3N0YXIsIGJldGEwID0gYmV0YTBfc3RhciwgYmV0YTEgPSB0aGV0YVtpLTEsIDNdKSAtCiAgICAgIGRub3JtKHggPSB0aGV0YVsoaS0xKSwzXSwgbWVhbiA9IGJldGExX3N0YXIsIHNkID0gMC4xLCBsb2cgPSBUUlVFKQojIyAhISEgaXMgdGhlIG1lYW4gY29ycmVjdD8gU2hvdWxkICB3ZSB1c2UgaW4gZGVuX2xvZ3IgbWVhbiA9IHNpZ21hU1Ffc3RhciBvciBtZWFuID0gdGhldGFbaS0xLCAxXT8gCiMgY2hlY2sgYy1jb21wb25lbnR3aXNlLW1oIHYxIGxpbmUgNjgsIDcyCiAgICBsb2dyID0gbnVtX2xvZ3IgLSBkZW5fbG9ncgogICAgaWYgKGxvZyhydW5pZigxKSkgPD0gbWluKGxvZ3IsIDApKSB7CiAgICAgIHRoZXRhW2ksM10gPSBiZXRhMV9zdGFyCiAgICB9IGVsc2UgewogICAgICB0aGV0YVtpLDNdID0gdGhldGFbKGkgLSAxKSwgM10KICAgIH0KICAgIAogIH0KICByZXR1cm4odGhldGEpCn0KCmBgYAoKYGBge3IgUHJvYmxlbSAxIG1oIFJ1bn0KCkIgPSAxMF41IAprZWVwID0gKEIvMiArIDEpOihCICsgMSkKY2hhaW4xID0gbWgoQiwgdGhldGFfc3RhcnQgPSBjKDAuMSwgLTEsIC0xKSkKY2hhaW4yID0gbWgoQiwgdGhldGFfc3RhcnQgPSBjKDAuMywgMCwgMCkpCmNoYWluMyA9IG1oKEIsIHRoZXRhX3N0YXJ0ID0gYygwLjUsIC0xLCAxKSkKY2hhaW40ID0gbWgoQiwgdGhldGFfc3RhcnQgPSBjKDAuMiwgMSwgMSkpCgptYyA9IG1jbWMubGlzdChtY21jKGNoYWluMVtrZWVwLF0pLCBtY21jKGNoYWluMltrZWVwLF0pLAogICAgICAgICAgICAgICAgIG1jbWMoY2hhaW4zW2tlZXAsXSksIG1jbWMoY2hhaW40W2tlZXAsXSkpCnN1bW1hcnkobWMpCgpgYGAKCmBgYHtyIFByb2JsZW0gMiB0byBhc3Nlc3MgdGhlIGNvbnZlcmdlbmNlfQoKa2VlcCA9IChCLzIgKyAxOTAwMSk6KEIvMiArIDIwMDAxKQoKbWMgPSBtY21jLmxpc3QobWNtYyhjaGFpbjFba2VlcCxdKSwgbWNtYyhjaGFpbjJba2VlcCxdKSwKICAgICAgICAgICAgICAgICBtY21jKGNoYWluM1trZWVwLF0pLCBtY21jKGNoYWluNFtrZWVwLF0pKQpjb2RhOjp0cmFjZXBsb3QobWMpCmNvZGE6OmF1dG9jb3JyLnBsb3QobWMsIGxhZy5tYXggPSAxMDAsIGF1dG8ubGF5b3V0ID0gVFJVRSkKZ2VsbWFuLmRpYWcobWMsIGF1dG9idXJuaW4gPSBGQUxTRSkKZ2V3ZWtlLmRpYWcobWMpCgpgYGAKCmBgYHtyIFByb2JsZW0gMyBNb2RlbCBBLCBpbmNsdWRlPUZBTFNFfQpzdGFuX21vZF9BID0gIgpkYXRhIHsKICBpbnQgVDsKICB2ZWN0b3JbVF0gbHByOyAgICAvLyBsb2cgcHJpY2UgCiAgdmVjdG9yW1RdIGxjYXI7ICAgLy8gbG9nIGNhcmF0IAogIHJlYWw8bG93ZXI9MD4gdjsgIC8vIHNhbXBsZSB2YXJpYW5jZSBvZiBsb2cgcHJpY2UgCgp9CgpwYXJhbWV0ZXJzIHsKICByZWFsPGxvd2VyPTA+IHNpZ21hU1E7CiAgcmVhbCBiZXRhMDsKICByZWFsIGJldGExOwp9Cgp0cmFuc2Zvcm1lZCBwYXJhbWV0ZXJzIHsKICB2ZWN0b3JbVF0gbXU7CiAgZm9yIChpIGluIDE6VCkgewogICAgbXVbaV0gPSBiZXRhMCArIGJldGExICogbGNhcltpXTsKICB9Cgp9Cgptb2RlbCB7CiAgc2lnbWFTUSB+IGludl9nYW1tYSgwLjAxLCAwLjAxKTsKICBiZXRhMCB+IG5vcm1hbCgwLCAxMDApOwogIGJldGExIH4gbm9ybWFsKDAsIDEwMCk7CgogIGZvciAoaSBpbiAxOlQpCiAgICBscHJbaV0gfiBub3JtYWwobXVbaV0sIHNxcnQoc2lnbWFTUSkpOwp9CmdlbmVyYXRlZCBxdWFudGl0aWVzIHsKICB2ZWN0b3JbVF0gbG9nX2xpazsKICB2ZWN0b3JbVF0geV9yZXA7CgogIC8vdmVjdG9yW1RdIGxpazsKICByZWFsIFJic3E7ICAgICAgICAgICAgICAvLyBnb29kbmVzcy1vZi1maXQKICBSYnNxID0gMSAtIHNpZ21hU1EvdjsKICAKICBmb3IgKGkgaW4gMTpUKSB7CiAgICBsb2dfbGlrW2ldID0gbm9ybWFsX2xwZGYobHByW2ldIHwgbXVbaV0sIHNxcnQoc2lnbWFTUSkpOwogICAgeV9yZXBbaV0gPSBub3JtYWxfcm5nKG11W2ldLCBzcXJ0KHNpZ21hU1EpKTsKICAgIC8vbGlrW2ldID0gZXhwKGxvZ19saWtbaV0pOwogIH0KCn0KIgoKVCA8LSBsZW5ndGgoZGYkbHByaWNlKQp2ID0gdmFyKGRmJGxwcmljZSkKc2V0LnNlZWQoOTApCgptb2RlbF9uYW1lID0gIk1vZGVsX0FfIgpzdGFuX2RhdF9BID0gbGlzdChUID0gVCwgbHByID0gZGYkbHByaWNlLCBsY2FyID0gZGYkbGNhcmF0LCB2ID0gdikKcl9maXRfQSA9IHN0YW4obW9kZWxfY29kZSA9IHN0YW5fbW9kX0EsIGRhdGEgPSBzdGFuX2RhdF9BLCAKICAgICAgICAgICAgIGl0ZXIgPSA1MDAwLCBjaGFpbnMgPSA0LCAKICAgICAgICAgICAgIGNvbnRyb2wgPSBsaXN0KG1heF90cmVlZGVwdGggPSAyMSkpCgojZmlsZV9uYW1lID0gcGFzdGUobW9kZWxfbmFtZSwgYXMuY2hhcmFjdGVyKFN5cy5EYXRlKCkpLCAiLnJkYSIsIHNlcD0iIikKI3NldHdkKCcvVXNlcnMvQU0vRG9jdW1lbnRzL19DVSBNYXN0ZXJzLzIwMjAgZmFsbCBCYXllc2lhbl83MzkzL2NvZGUvQmF5ZXNpYW5fU3RhdGlzdGljc19DbGFzc19Db2RlL0V4YW1fY29kZScpCiNyZWFkUkRTKHJfZml0LCBmaWxlID0gIiIpIAojc2F2ZVJEUyhyX2ZpdCwgZmlsZSA9IGZpbGVfbmFtZSwgY29tcHJlc3MgPSAieHoiKSAKCnN1bW1hcnkocl9maXRfQSwgcGFycyA9IGMoInNpZ21hU1EiLCAiYmV0YTAiLCAiYmV0YTEiKSwgcHJvYiA9IGMoMC4wMjUsIDAuOTc1KSkkc3VtbWFyeQoKI3NzbyA8LSBsYXVuY2hfc2hpbnlzdGFuKHJfZml0X0EpCgpgYGAKCmBgYHtyIFByb2JsZW0gMyBNb2RlbCBCLCBpbmNsdWRlPUZBTFNFfQpzdGFuX21vZF9CID0gIgpkYXRhIHsKICBpbnQgVDsKICB2ZWN0b3JbVF0gcHI7ICAgICAvLyBwcmljZSAKICB2ZWN0b3JbVF0gY2FyOyAgICAvLyBjYXJhdCAKICByZWFsPGxvd2VyPTA+IHY7ICAvLyBzYW1wbGUgdmFyaWFuY2Ugb2YgbG9nIHByaWNlIAoKfQoKcGFyYW1ldGVycyB7CiAgcmVhbDxsb3dlcj0wPiBzaWdtYVNROwogIHJlYWwgYmV0YTA7CiAgcmVhbCBiZXRhMTsKICByZWFsIGJldGEyOwp9Cgp0cmFuc2Zvcm1lZCBwYXJhbWV0ZXJzIHsKICB2ZWN0b3JbVF0gbXU7CiAgZm9yIChpIGluIDE6VCkgewogICAgbXVbaV0gPSBiZXRhMCArIGJldGExICogY2FyW2ldICsgYmV0YTIgKiBjYXJbaV0gKiBjYXJbaV07CiAgfQoKfQoKbW9kZWwgewogIHNpZ21hU1EgfiBpbnZfZ2FtbWEoMC4wMSwgMC4wMSk7CiAgYmV0YTAgfiBub3JtYWwoMCwgMTAwKTsKICBiZXRhMSB+IG5vcm1hbCgwLCAxMDApOwogIGJldGEyIH4gbm9ybWFsKDAsIDEwMCk7CgogIGZvciAoaSBpbiAxOlQpCiAgICBwcltpXSB+IG5vcm1hbChtdVtpXSwgc3FydChzaWdtYVNRKSk7Cn0KCmdlbmVyYXRlZCBxdWFudGl0aWVzIHsKICB2ZWN0b3JbVF0geV9yZXA7CiAgdmVjdG9yW1RdIGxvZ19saWs7CiAgLy92ZWN0b3JbVF0gbGlrOwogIHJlYWwgUmJzcTsgICAgICAgICAgICAgIC8vIGdvb2RuZXNzLW9mLWZpdAogIFJic3EgPSAxIC0gc2lnbWFTUS92OwogIAogIGZvciAoaSBpbiAxOlQpIHsKICAgIGxvZ19saWtbaV0gPSBub3JtYWxfbHBkZihwcltpXSB8IG11W2ldLCBzcXJ0KHNpZ21hU1EpKTsKICAgIHlfcmVwW2ldID0gbm9ybWFsX3JuZyhtdVtpXSwgc3FydChzaWdtYVNRKSk7CgogICAgLy9saWtbaV0gPSBleHAobG9nX2xpa1tpXSk7CiAgfQoKfQoiCgpUIDwtIGxlbmd0aChkZiRwcmljZSkKdiA9IHZhcihkZiRwcmljZSkKc2V0LnNlZWQoOTEpCgptb2RlbF9uYW1lID0gIk1vZGVsX0JfIgpzdGFuX2RhdF9CID0gbGlzdChUID0gVCwgcHIgPSBkZiRwcmljZSwgY2FyID0gZGYkY2FyYXQsIHYgPSB2KQpyX2ZpdF9CID0gc3Rhbihtb2RlbF9jb2RlID0gc3Rhbl9tb2RfQiwgZGF0YSA9IHN0YW5fZGF0X0IsIAogICAgICAgICAgICAgaXRlciA9IDUwMDAsIGNoYWlucyA9IDQsIAogICAgICAgICAgICAgY29udHJvbCA9IGxpc3QobWF4X3RyZWVkZXB0aCA9IDIxKSkKCiNmaWxlX25hbWUgPSBwYXN0ZShtb2RlbF9uYW1lLCBhcy5jaGFyYWN0ZXIoU3lzLnRpbWUoKSksICIucmRhIiwgc2VwPSIiKQojc2V0d2QoJy9Vc2Vycy9BTS9Eb2N1bWVudHMvX0NVIE1hc3RlcnMvMjAyMCBmYWxsIEJheWVzaWFuXzczOTMvY29kZS9CYXllc2lhbl9TdGF0aXN0aWNzX0NsYXNzX0NvZGUvRXhhbV9jb2RlJykKI3NhdmVSRFMocl9maXQsIGZpbGUgPSBmaWxlX25hbWUsIGNvbXByZXNzID0gInh6IikgCnN1bW1hcnkocl9maXRfQiwgcGFycyA9IGMoInNpZ21hU1EiLCAiYmV0YTAiLCAiYmV0YTEiLCAiYmV0YTIiKSwgcHJvYiA9IGMoMC4wMjUsIDAuOTc1KSkkc3VtbWFyeQoKc3NvIDwtIGxhdW5jaF9zaGlueXN0YW4ocl9maXRfQikKCmBgYAoKYGBge3IgUHJvYmxlbSA0IHBwX2NoZWNrIE1vZGVsIEF9CiMgZXh0cmFjdGluZyB0aGUgTUNNQyBzYW1wbGVzIGZyb20gdGhlIHNvZGFfZml0IG9iamVjdApzYW1wbGVzID0gZXh0cmFjdChyX2ZpdF9BKQpuY3ljbGVzID0gbGVuZ3RoKHNhbXBsZXNbWzFdXSkKVCA8LSBsZW5ndGgoZGYkcHJpY2UpCgojIGVhY2ggcm93IG9mIHlyZXAgaXMgYSBzYW1wbGUgZnJvbSB0aGUgcHAgZGlzdHJpYnV0aW9uCnlyZXAgPSBtYXRyaXgoMCwgbmNvbCA9IFQsIG5yb3cgPSBuY3ljbGVzKQpmb3IgKGkgaW4gc2VxX2xlbihUKSkgewogIG11aSA9IHNhbXBsZXMkYmV0YTAgKyBzYW1wbGVzJGJldGExICogZGYkbGNhcmF0W2ldCiAgeXJlcFssIGldID0gcm5vcm0obmN5Y2xlcywgbWVhbiA9IG11aSwgc2QgPSBzcXJ0KHNhbXBsZXMkc2lnbWFTUSkpCn0KCiMgYXBwcm94aW1hdGUgcG9zdGVyaW9yIHByZWRpY3RpdmUgZGVuc2l0eSBmb3IgYW4gb2JzZXJ2YXRpb24KIyB3aXRoIHRoZSBzYW1lIGNvdmFyaWF0ZXMgYXMgb2JzZXJ2YXRpb24gMQojcGxvdChkZW5zaXR5KHlyZXBbLDFdKSkKY29sb3Jfc2NoZW1lX3NldCgicmVkIikKCiMgcG9zdGVyaW9yIHByZWRpY3RpdmUgY2hlY2sKeSA9IGRmJGxwcmljZQpwcGNfaGlzdCh5LCB5cmVwW3NhbXBsZSgxOm5jeWNsZXMsIDgpLF0pICsgZ2d0aXRsZSgiTW9kZWwgQTogaGlzdG9ncmFtIGNvbXBhcmluZyB0aGUgcmVzcG9uc2UgdG8gOCByZXBsaWNhdGVkIGRhdGEgc2V0cyIpCgoKIyBzY2F0dGVycGxvdCBvZiB5IHZzIGF2ZXJhZ2UgeXJlcApwcGNfc2NhdHRlcl9hdmcoeSwgeXJlcCkgKyBnZ3RpdGxlKCJNb2RlbCBBOiBzY2F0dGVycGxvdCBjb21wYXJpbmcgdGhlIHJlc3BvbnNlIHRvIHRoZSBhdmVyYWdlIHJlcGxpY2F0ZWQgZGF0YSBzZXQiKQoKcHBjX2RlbnNfb3ZlcmxheSh5LCB5cmVwID0geXJlcFtzYW1wbGUoMTpuY3ljbGVzLCAxMDApLF0pICsgZ2d0aXRsZSgiTW9kZWwgQTogZGVuc2l0eSBwbG90IGNvbXBhcmluZyB0aGUgZGVuc2l0eSBvZiB0aGUgcmVzcG9uc2UgdG8gdGhlIGRlbnNpdGllcyBvZiAxMDAgcmVwbGljYXRlZCBkYXRhIHNldHMiKQoKcHBjX2Vycm9yX3NjYXR0ZXJfYXZnKHksIHlyZXAgPSB5cmVwKSArIGdndGl0bGUoIk1vZGVsIEE6IHNjYXR0ZXIgcGxvdCBjb21wYXJpbmcgdGhlIHJlc3BvbnNlIHRvIHRoZSBhdmVyYWdlIHByZWRpY3RpdmUgZXJyb3IiKQoKcHBjX2ludGVydmFscyh5ID0geSwgeXJlcCA9IHlyZXAsIHggPSBkZiRsY2FyYXQpICsgZ2d0aXRsZSgiTW9kZWwgQTogUG9zdGVyaW9yIHByZWRpY3RpdmUgaW50ZXJ2YWxzIHRvIHRoZSBvYnNlcnZlZCBkYXRhIHZhbHVlcyIpICsgeGxhYigibG9nIGNhcmF0IikgKyB5bGFiKCJsb2cgcHJpY2UiKQpgYGAKCmBgYHtyIFByb2JsZW0gNCBwcF9jaGVjayBNb2RlbCBCfQojIGV4dHJhY3RpbmcgdGhlIE1DTUMgc2FtcGxlcyBmcm9tIHRoZSBzb2RhX2ZpdCBvYmplY3QKc2FtcGxlcyA9IGV4dHJhY3Qocl9maXRfQikKbmN5Y2xlcyA9IGxlbmd0aChzYW1wbGVzW1sxXV0pClQgPC0gbGVuZ3RoKGRmJHByaWNlKQoKIyBlYWNoIHJvdyBvZiB5cmVwIGlzIGEgc2FtcGxlIGZyb20gdGhlIHBwIGRpc3RyaWJ1dGlvbgp5cmVwID0gbWF0cml4KDAsIG5jb2wgPSBULCBucm93ID0gbmN5Y2xlcykKZm9yIChpIGluIHNlcV9sZW4oVCkpIHsKICBtdWkgPSBzYW1wbGVzJGJldGEwICsgc2FtcGxlcyRiZXRhMSAqIGRmJGNhcmF0W2ldICsgc2FtcGxlcyRiZXRhMiAqIGRmJGNhcmF0W2ldICogZGYkY2FyYXRbaV0KICB5cmVwWywgaV0gPSBybm9ybShuY3ljbGVzLCBtZWFuID0gbXVpLCBzZCA9IHNxcnQoc2FtcGxlcyRzaWdtYVNRKSkKfQoKIyBhcHByb3hpbWF0ZSBwb3N0ZXJpb3IgcHJlZGljdGl2ZSBkZW5zaXR5IGZvciBhbiBvYnNlcnZhdGlvbgojIHdpdGggdGhlIHNhbWUgY292YXJpYXRlcyBhcyBvYnNlcnZhdGlvbiAxCiNwbG90KGRlbnNpdHkoeXJlcFssMV0pKQpjb2xvcl9zY2hlbWVfc2V0KCJncmVlbiIpCgojIHBvc3RlcmlvciBwcmVkaWN0aXZlIGNoZWNrCnkgPSBkZiRwcmljZQpwcGNfaGlzdCh5LCB5cmVwW3NhbXBsZSgxOm5jeWNsZXMsIDgpLF0pICsgZ2d0aXRsZSgiTW9kZWwgQjogaGlzdG9ncmFtIGNvbXBhcmluZyB0aGUgcmVzcG9uc2UgdG8gOCByZXBsaWNhdGVkIGRhdGEgc2V0cyIpCgoKIyBzY2F0dGVycGxvdCBvZiB5IHZzIGF2ZXJhZ2UgeXJlcApwcGNfc2NhdHRlcl9hdmcoeSwgeXJlcCkgKyBnZ3RpdGxlKCJNb2RlbCBCOiBzY2F0dGVycGxvdCBjb21wYXJpbmcgdGhlIHJlc3BvbnNlIHRvIHRoZSBhdmVyYWdlIHJlcGxpY2F0ZWQgZGF0YSBzZXQiKQoKcHBjX2RlbnNfb3ZlcmxheSh5LCB5cmVwID0geXJlcFtzYW1wbGUoMTpuY3ljbGVzLCAxMDApLF0pICsgZ2d0aXRsZSgiTW9kZWwgQjogZGVuc2l0eSBwbG90IGNvbXBhcmluZyB0aGUgZGVuc2l0eSBvZiB0aGUgcmVzcG9uc2UgdG8gdGhlIGRlbnNpdGllcyBvZiAxMDAgcmVwbGljYXRlZCBkYXRhIHNldHMiKQoKcHBjX2Vycm9yX3NjYXR0ZXJfYXZnKHksIHlyZXAgPSB5cmVwKSArIGdndGl0bGUoIk1vZGVsIEI6IHNjYXR0ZXIgcGxvdCBjb21wYXJpbmcgdGhlIHJlc3BvbnNlIHRvIHRoZSBhdmVyYWdlIHByZWRpY3RpdmUgZXJyb3IiKQoKcHBjX2ludGVydmFscyh5ID0geSwgeXJlcCA9IHlyZXAsIHggPSBkZiRjYXJhdCkgKyBnZ3RpdGxlKCJNb2RlbCBCOiBQb3N0ZXJpb3IgcHJlZGljdGl2ZSBpbnRlcnZhbHMgdG8gdGhlIG9ic2VydmVkIGRhdGEgdmFsdWVzIikgKyB4bGFiKCJjYXJhdCIpICsgeWxhYigicHJpY2UiKQpgYGAKCmBgYHtyIFByb2JsZW0gNSBNb2RlbCBDLCBpbmNsdWRlPUZBTFNFfQpzdGFuX21vZF9DID0gIgpkYXRhIHsKICBpbnQgVDsKICB2ZWN0b3JbVF0gcHI7ICAgIC8vICBwcmljZSAKICB2ZWN0b3JbVF0gbGNhcjsgICAvLyBsb2cgY2FyYXQgCiAgcmVhbDxsb3dlcj0wPiB2OyAgLy8gc2FtcGxlIHZhcmlhbmNlIG9mIGxvZyBwcmljZSAKCn0KCnBhcmFtZXRlcnMgewogIHJlYWw8bG93ZXI9MD4gc2lnbWFTUTsKICByZWFsIGJldGEwOwogIHJlYWwgYmV0YTE7Cn0KCnRyYW5zZm9ybWVkIHBhcmFtZXRlcnMgewogIHZlY3RvcltUXSBtdTsKICBmb3IgKGkgaW4gMTpUKSB7CiAgICBtdVtpXSA9IGJldGEwICsgYmV0YTEgKiBsY2FyW2ldOwogIH0KCn0KCm1vZGVsIHsKICBzaWdtYVNRIH4gaW52X2dhbW1hKDAuMDEsIDAuMDEpOwogIGJldGEwIH4gbm9ybWFsKDAsIDEwMCk7CiAgYmV0YTEgfiBub3JtYWwoMCwgMTAwKTsKCiAgZm9yIChpIGluIDE6VCkKICAgIHByW2ldIH4gbG9nbm9ybWFsKG11W2ldLCBzcXJ0KHNpZ21hU1EpKTsKfQpnZW5lcmF0ZWQgcXVhbnRpdGllcyB7CiAgdmVjdG9yW1RdIGxvZ19saWs7CiAgdmVjdG9yW1RdIHlfcmVwOwoKICAvL3ZlY3RvcltUXSBsaWs7CiAgcmVhbCBSYnNxOyAgICAgICAgICAgICAgLy8gZ29vZG5lc3Mtb2YtZml0CiAgUmJzcSA9IDEgLSBzaWdtYVNRL3Y7CiAgCiAgZm9yIChpIGluIDE6VCkgewogICAgbG9nX2xpa1tpXSA9IGxvZ25vcm1hbF9scGRmKHByW2ldIHwgbXVbaV0sIHNxcnQoc2lnbWFTUSkpOwogICAgeV9yZXBbaV0gPSBsb2dub3JtYWxfcm5nKG11W2ldLCBzcXJ0KHNpZ21hU1EpKTsKICAgIC8vbGlrW2ldID0gZXhwKGxvZ19saWtbaV0pOwogIH0KCn0KIgoKVCA8LSBsZW5ndGgoZGYkcHJpY2UpCnYgPSB2YXIoZGYkcHJpY2UpCnNldC5zZWVkKDkyKQoKbW9kZWxfbmFtZSA9ICJNb2RlbF9DXyIKc3Rhbl9kYXRfQyA9IGxpc3QoVCA9IFQsIHByID0gZGYkcHJpY2UsIGxjYXIgPSBkZiRsY2FyYXQsIHYgPSB2KQpyX2ZpdF9DID0gc3Rhbihtb2RlbF9jb2RlID0gc3Rhbl9tb2RfQywgZGF0YSA9IHN0YW5fZGF0X0MsIAogICAgICAgICAgICAgaXRlciA9IDUwMDAsIGNoYWlucyA9IDQsIAogICAgICAgICAgICAgY29udHJvbCA9IGxpc3QobWF4X3RyZWVkZXB0aCA9IDIxKSkKCiNmaWxlX25hbWUgPSBwYXN0ZShtb2RlbF9uYW1lLCBhcy5jaGFyYWN0ZXIoU3lzLkRhdGUoKSksICIucmRhIiwgc2VwPSIiKQojc2V0d2QoJy9Vc2Vycy9BTS9Eb2N1bWVudHMvX0NVIE1hc3RlcnMvMjAyMCBmYWxsIEJheWVzaWFuXzczOTMvY29kZS9CYXllc2lhbl9TdGF0aXN0aWNzX0NsYXNzX0NvZGUvRXhhbV9jb2RlJykKI3JlYWRSRFMocl9maXQsIGZpbGUgPSAiIikgCiNzYXZlUkRTKHJfZml0LCBmaWxlID0gZmlsZV9uYW1lLCBjb21wcmVzcyA9ICJ4eiIpIAoKc3VtbWFyeShyX2ZpdF9DLCBwYXJzID0gYygic2lnbWFTUSIsICJiZXRhMCIsICJiZXRhMSIpLCBwcm9iID0gYygwLjAyNSwgMC45NzUpKSRzdW1tYXJ5Cgojc3NvIDwtIGxhdW5jaF9zaGlueXN0YW4ocl9maXRfQSkKCmBgYAoKYGBge3IgUHJvYmxlbSA1IE1vZGVsIEQsIGluY2x1ZGU9RkFMU0V9CnN0YW5fbW9kX0QgPSAiCmRhdGEgewogIGludCBUOwogIHZlY3RvcltUXSBwcjsgICAgIC8vICBwcmljZSAKICB2ZWN0b3JbVF0gbGNhcjsgICAvLyBsb2cgY2FyYXQgCiAgdmVjdG9yW1RdIG5vdGljOyAgLy8gIG5vdGljZWFibGUgCgogIHJlYWw8bG93ZXI9MD4gdjsgIC8vIHNhbXBsZSB2YXJpYW5jZSBvZiBsb2cgcHJpY2UgCgp9CgpwYXJhbWV0ZXJzIHsKICByZWFsPGxvd2VyPTA+IHNpZ21hU1E7CiAgcmVhbCBiZXRhMDsKICByZWFsIGJldGExOwogIHJlYWwgYmV0YTI7Cn0KCnRyYW5zZm9ybWVkIHBhcmFtZXRlcnMgewogIHZlY3RvcltUXSBtdTsKICBmb3IgKGkgaW4gMTpUKSB7CiAgICBtdVtpXSA9IGJldGEwICsgYmV0YTEgKiBsY2FyW2ldICsgYmV0YTIgKiBub3RpY1tpXTsKICB9Cgp9Cgptb2RlbCB7CiAgc2lnbWFTUSB+IGludl9nYW1tYSgwLjAxLCAwLjAxKTsKICBiZXRhMCB+IG5vcm1hbCgwLCAxMDApOwogIGJldGExIH4gbm9ybWFsKDAsIDEwMCk7CiAgYmV0YTIgfiBub3JtYWwoMCwgMTAwKTsKCgogIGZvciAoaSBpbiAxOlQpCiAgICBwcltpXSB+IGxvZ25vcm1hbChtdVtpXSwgc3FydChzaWdtYVNRKSk7Cn0KZ2VuZXJhdGVkIHF1YW50aXRpZXMgewogIHZlY3RvcltUXSBsb2dfbGlrOwogIHZlY3RvcltUXSB5X3JlcDsKCiAgLy92ZWN0b3JbVF0gbGlrOwogIHJlYWwgUmJzcTsgICAgICAgICAgICAgIC8vIGdvb2RuZXNzLW9mLWZpdAogIFJic3EgPSAxIC0gc2lnbWFTUS92OwogIAogIGZvciAoaSBpbiAxOlQpIHsKICAgIGxvZ19saWtbaV0gPSBsb2dub3JtYWxfbHBkZihwcltpXSB8IG11W2ldLCBzcXJ0KHNpZ21hU1EpKTsKICAgIHlfcmVwW2ldID0gbG9nbm9ybWFsX3JuZyhtdVtpXSwgc3FydChzaWdtYVNRKSk7CiAgICAvL2xpa1tpXSA9IGV4cChsb2dfbGlrW2ldKTsKICB9Cgp9CiIKClQgPC0gbGVuZ3RoKGRmJHByaWNlKQp2ID0gdmFyKGRmJHByaWNlKQpzZXQuc2VlZCg5MykKCm1vZGVsX25hbWUgPSAiTW9kZWxfRF8iCnN0YW5fZGF0X0QgPSBsaXN0KFQgPSBULCBwciA9IGRmJHByaWNlLCBsY2FyID0gZGYkbGNhcmF0LCAKICAgICAgICAgICAgICAgICAgbm90aWMgPSBhcy5pbnRlZ2VyKGRmJGluY2x1c2lvbnMgPT0gIm5vdGljZWFibGUiKSAsIHYgPSB2KQpyX2ZpdF9EID0gc3Rhbihtb2RlbF9jb2RlID0gc3Rhbl9tb2RfRCwgZGF0YSA9IHN0YW5fZGF0X0QsIAogICAgICAgICAgICAgaXRlciA9IDUwMDAsIGNoYWlucyA9IDQsIAogICAgICAgICAgICAgY29udHJvbCA9IGxpc3QobWF4X3RyZWVkZXB0aCA9IDIxKSkKCiNmaWxlX25hbWUgPSBwYXN0ZShtb2RlbF9uYW1lLCBhcy5jaGFyYWN0ZXIoU3lzLkRhdGUoKSksICIucmRhIiwgc2VwPSIiKQojc2V0d2QoJy9Vc2Vycy9BTS9Eb2N1bWVudHMvX0NVIE1hc3RlcnMvMjAyMCBmYWxsIEJheWVzaWFuXzczOTMvY29kZS9CYXllc2lhbl9TdGF0aXN0aWNzX0NsYXNzX0NvZGUvRXhhbV9jb2RlJykKI3JlYWRSRFMocl9maXQsIGZpbGUgPSAiIikgCiNzYXZlUkRTKHJfZml0LCBmaWxlID0gZmlsZV9uYW1lLCBjb21wcmVzcyA9ICJ4eiIpIAoKc3VtbWFyeShyX2ZpdF9ELCBwYXJzID0gYygic2lnbWFTUSIsICJiZXRhMCIsICJiZXRhMSIsICJiZXRhMiIpLCBwcm9iID0gYygwLjAyNSwgMC45NzUpKSRzdW1tYXJ5Cgojc3NvIDwtIGxhdW5jaF9zaGlueXN0YW4ocl9maXRfQSkKCmBgYAoKYGBge3IgUHJvYmxlbSA2YSBXQUlDL0xPT0lDLCBpbmNsdWRlPUZBTFNFfQoKd2FpY19sb28gPSBmdW5jdGlvbiAoZml0KSB7CiAgbG9nX2xpayA9IGV4dHJhY3RfbG9nX2xpayhmaXQsIG1lcmdlX2NoYWlucyA9IEZBTFNFKQogIHJfZWZmID0gZXhwKHJlbGF0aXZlX2VmZihsb2dfbGlrKSkKICB3bCA9IGMod2FpYyhsb2dfbGlrKSR3YWljLCBsb28obG9nX2xpaywgcl9lZmYgPSByX2VmZikkbG9vaWMpCiAgcmV0dXJuKHdsKQp9CndsID0gYXJyYXkoMCwgZGltPWMoMiwzKSwgCiAgICAgICAgICAgZGltbmFtZXMgPSBsaXN0KGMoIldBSUMiLCAiTE9PSUMiKSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgIGMoIk1vZGVsIEIiLCJNb2RlbCBDIiwiTW9kZWwgRCIpKSkKd2xbMSwxXSA9IHdhaWNfbG9vKHJfZml0X0IpWzFdCndsWzIsMV0gPSB3YWljX2xvbyhyX2ZpdF9CKVsyXQp3bFsxLDJdID0gd2FpY19sb28ocl9maXRfQylbMV0Kd2xbMiwyXSA9IHdhaWNfbG9vKHJfZml0X0MpWzJdCndsWzEsM10gPSB3YWljX2xvbyhyX2ZpdF9EKVsxXQp3bFsyLDNdID0gd2FpY19sb28ocl9maXRfRClbMl0KCndsCmBgYAoKYGBge3IgUHJvYmxlbSA2YiBlZmZlY3RzIHBsb3RzfQpzYW1wbGVzID0gZXh0cmFjdChyX2ZpdF9EKQpuY3ljbGVzID0gbGVuZ3RoKHNhbXBsZXNbWzFdXSkKVCA8LSBsZW5ndGgoZGYkcHJpY2UpCmRmJG5vdGljID0gYXMuaW50ZWdlcihkZiRpbmNsdXNpb25zID09ICJub3RpY2VhYmxlIikKCiMgZWFjaCByb3cgb2YgeXJlcCBpcyBhIHNhbXBsZSBmcm9tIHRoZSBwcCBkaXN0cmlidXRpb24KeXJlcCA9IG1hdHJpeCgwLCBuY29sID0gVCwgbnJvdyA9IG5jeWNsZXMpCmZvciAoaSBpbiBzZXFfbGVuKFQpKSB7CiAgbXVpID0gc2FtcGxlcyRiZXRhMCArIHNhbXBsZXMkYmV0YTEgKiBkZiRsY2FyYXRbaV0gKyBzYW1wbGVzJGJldGEyICogZGYkbm90aWNbaV0KICB5cmVwWywgaV0gPSBybG5vcm0obmN5Y2xlcywgbWVhbiA9IG11aSwgc2QgPSBzcXJ0KHNhbXBsZXMkc2lnbWFTUSkpCn0KCmNvbG9yX3NjaGVtZV9zZXQoImJsdWUiKQpwcGNfaW50ZXJ2YWxzX2dyb3VwZWQoeSA9IGRmJHByaWNlLCB5cmVwID0geXJlcCwgeCA9IGRmJGNhcmF0LCBwcm9iID0gMC41LCBwcm9iX291dGVyID0gLjk1LCBncm91cCA9IGRmJG5vdGljKSArIGdndGl0bGUoIk1vZGVsIEQ6IFBvc3RlcmlvciBwcmVkaWN0aXZlIGludGVydmFscyB0byB0aGUgb2JzZXJ2ZWQgZGF0YSB2YWx1ZXMiKSArIHhsYWIoImxjYXJhdCIpICsgeWxhYigicHJpY2UiKQoKCmZpbmRfZXAgPSBmdW5jdGlvbih5cmVwLCBtaW4gPSAwLjEsIG1heCA9IDIuNCwgc3RlcD0wLjEpIHsKICBlcCA9IGFycmF5KDAsIGRpbT1jKChtYXgtbWluK3N0ZXApL3N0ZXAsIDUpLCBkaW1uYW1lcyA9IGxpc3QoYygpLCBjKCJjYXJhdF9taW4iLCAiY2FyYXRfbWF4IiwgIjIuNSUiLCJNZWFuIiwiOTcuNSUiKSkpCiAgbj0wCiAgZm9yIChpIGluIHNlcShtaW4rc3RlcCwgbWF4LCBieT1zdGVwKSkgewogICAgbj1uKzEKICAgIHN1YnNldCA9ICgoZGYkY2FyYXQ+KGktc3RlcCkpICYgKGRmJGNhcmF0PD1pKSkKICAgIHlyZXBfc3Vic2V0ID0geXJlcFssc3Vic2V0XQogICAgZXBbbiwxXSA9IGktc3RlcAogICAgZXBbbiwyXSA9IGkKICAgIHEgPSBxdWFudGlsZSh5cmVwX3N1YnNldCwgcHJvYnMgPSBjKDAuMDI1LCAwLjk3NSkpCiAgICBlcFtuLDNdID0gcVtbMV1dCiAgICBlcFtuLDRdID0gbWVhbih5cmVwX3N1YnNldCkKICAgIGVwW24sNV0gPSBxW1syXV0KICB9CnJldHVybihlcCkKfQoKZmluZF9lcCh5cmVwKQp0ZW1wID0gcHBjX2ludGVydmFsc19kYXRhKHkgPSBkZiRwcmljZSwgeXJlcCA9IHlyZXAsIHggPSBkZiRsY2FyYXQsIHByb2IgPSAwLjUsIHByb2Jfb3V0ZXIgPSAuOTUsIGdyb3VwID0gZGYkbm90aWMpCgpgYGAKCmBgYHtyIFByb2JsZW0gN19hIE1vZGVsIEUsIGluY2x1ZGU9RkFMU0V9CgpzdGFuX21vZF9FID0gIgpkYXRhIHsKICBpbnQgVDsKICBpbnQgZGVhdGhzW1RdOyAvLyBDT1ZJRC0xOSBkZWF0aHMKICB2ZWN0b3JbVF0gbG9nX3BvcDsgICAgLy8gbG9nIG9mIFUuUy4gc3RhdGVzIHBvcHVsYXRpb24sCiAgdmVjdG9yW1RdIGluY29tZTsgLy8gbWVkaWFuIGluY29tZSAoVVNEKQogIHZlY3RvcltUXSBiczsgICAgIC8vIHBlcmNlbnRhZ2Ugb2YgdGhlIHBvcHVsYXRpb24gd2l0aCBiYWNoZWxvcuKAmXMgZGVncmVlcy4KfQoKcGFyYW1ldGVycyB7CiAgcmVhbCBiZXRhMDsKICByZWFsIGJldGExOwogIHJlYWwgYmV0YTI7Cn0KdHJhbnNmb3JtZWQgcGFyYW1ldGVycyB7CiAgcmVhbCBsb2dfbGFtYmRhW1RdOwogIGZvciAoaSBpbiAxOlQpIHsKICAgbG9nX2xhbWJkYVtpXSA9IGxvZ19wb3BbaV0gKyBiZXRhMCArIGJldGExICogaW5jb21lW2ldICsgYmV0YTIgKiBic1tpXTsKICB9Cgp9Cgptb2RlbCB7CiAgYmV0YTAgfiBub3JtYWwoMCwgNSk7CiAgYmV0YTEgfiBub3JtYWwoMCwgNSk7CiAgYmV0YTIgfiBub3JtYWwoMCwgNSk7CiAKICBmb3IgKGkgaW4gMTpUKSB7CiAgICBkZWF0aHNbaV0gfiBwb2lzc29uX2xvZyhsb2dfbGFtYmRhW2ldKTsKICB9Cn0KZ2VuZXJhdGVkIHF1YW50aXRpZXMgewogIHZlY3RvcltUXSBsb2dfbGlrOwoKICBmb3IgKGkgaW4gMTpUKSB7CiAgICBsb2dfbGlrW2ldID0gcG9pc3Nvbl9scG1mKGRlYXRoc1tpXSB8IGV4cChsb2dfbGFtYmRhW2ldKSk7CiAgfSAgCiAgCn0KCiIKClQgPC0gbGVuZ3RoKGRmX2NvdmlkJGJzKQpzZXQuc2VlZCg5NCkKCm1vZGVsX25hbWUgPSAiTW9kZWxfRV8iCnN0YW5fZGF0X0UgPSBsaXN0KFQgPSBULCBkZWF0aHMgPSBkZl9jb3ZpZCRkZWF0aHMsIGxvZ19wb3AgPSBsb2coZGZfY292aWQkcG9wdWxhdGlvbiksIAogICAgICAgICAgICAgICAgICBpbmNvbWUgPSBkZl9jb3ZpZCRpbmNvbWUsIGJzID0gZGZfY292aWQkYnMpCnJfZml0X0UgPSBzdGFuKG1vZGVsX2NvZGUgPSBzdGFuX21vZF9FLCBkYXRhID0gc3Rhbl9kYXRfRSwgCiAgICAgICAgICAgICBpdGVyID0gNTAwMDAsIGNoYWlucyA9IDIsIAogICAgICAgICAgICAgY29udHJvbCA9IGxpc3QobWF4X3RyZWVkZXB0aCA9IDIxKSkKCnN1bW1hcnkocl9maXRfRSwgcGFycyA9IGMoImJldGEwIiwgImJldGExIiwgImJldGEyIiksIHByb2IgPSBjKDAuMDI1LCAwLjk3NSkpJHN1bW1hcnkKCiNzc28gPC0gbGF1bmNoX3NoaW55c3RhbihyX2ZpdF9FKQpmaWxlX25hbWUgPSBwYXN0ZShtb2RlbF9uYW1lLCBhcy5jaGFyYWN0ZXIoU3lzLkRhdGUoKSksICIucmRhIiwgc2VwPSIiKQpzZXR3ZCgnL1VzZXJzL0FNL0RvY3VtZW50cy9fQ1UgTWFzdGVycy8yMDIwIGZhbGwgQmF5ZXNpYW5fNzM5My9jb2RlL0JheWVzaWFuX1N0YXRpc3RpY3NfQ2xhc3NfQ29kZS9FeGFtX2NvZGUnKQojcmVhZFJEUyhyX2ZpdCwgZmlsZSA9ICIiKSAKI3NhdmVSRFMocl9maXRfRSwgZmlsZSA9IGZpbGVfbmFtZSwgY29tcHJlc3MgPSAieHoiKSAKYGBgCgpgYGB7ciBQcm9ibGVtIDdfYiByc3RhbmFybSBNb2RlbCBFLCBpbmNsdWRlPUZBTFNFfQpkZl9jb3ZpZF9yc3RhbiA9IGRhdGEuZnJhbWUoaW50ID0gMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBkZWF0aHMgPSBkZl9jb3ZpZCRkZWF0aHMsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBsb2dfcG9wID0gbG9nKGRmX2NvdmlkJHBvcHVsYXRpb24pLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGluY29tZSA9IGRmX2NvdmlkJGluY29tZSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGJzID0gZGZfY292aWQkYnMpCmZpdF9yc3RhbmFybV9FID0gc3Rhbl9nbG0oZGVhdGhzIH4gaW50ICsgaW5jb21lICsgYnMsIAogICAgICAgICAgICAgICAgIG9mZnNldCA9IGxvZ19wb3AsCiAgICAgICAgICAgICAgICAgZmFtaWx5ID0gcG9pc3NvbihsaW5rID0gImxvZyIpLAogICAgICAgICAgICAgICAgIGRhdGEgPSBkZl9jb3ZpZF9yc3RhbiwKICAgICAgICAgICAgICAgICBpdGVyID0gMTAwMDAsIGNoYWlucyA9IDIpCgpzdW1tYXJ5KGZpdF9yc3RhbmFybV9FLCBkaWdpdHMgPSA2KQoKLTAuMDAwMDU0ICogMzA2OTcKKDI0MTkxIC0gMzk2ODIpICogLTAuMDAwMDU0CgooMTcuNTAgLSAzOS4wMCkgKiAwLjA2MDM0OApgYGAKCmBgYHtyIFByb2JsZW0gN19jIFNUQU4gTW9kZWwgRSBXQUlDIExPT0lDfQpjYXQoIlNUQU4gTW9kZWwgRSBXQUlDIiwgd2FpY19sb28ocl9maXRfRSlbMV0pCmNhdCgiU1RBTiBNb2RlbCBFIExPT0lDIiwgd2FpY19sb28ocl9maXRfRSlbMl0pCmBgYAoKYGBge3IgUHJvYmxlbSA3X2QgU1RBTiBNb2RlbCBFIHBwY30Kc2FtcGxlcyA9IGV4dHJhY3Qocl9maXRfRSkKbmN5Y2xlcyA9IGxlbmd0aChzYW1wbGVzW1sxXV0pClQgPC0gbGVuZ3RoKGRmX2NvdmlkJGJzKQoKIyBlYWNoIHJvdyBvZiB5cmVwIGlzIGEgc2FtcGxlIGZyb20gdGhlIHBwIGRpc3RyaWJ1dGlvbgp5cmVwID0gbWF0cml4KDAsIG5jb2wgPSBULCBucm93ID0gbmN5Y2xlcykKZm9yIChpIGluIHNlcV9sZW4oVCkpIHsKICBsb2dfbCA9IGxvZyhkZl9jb3ZpZCRwb3B1bGF0aW9uW2ldKSArCiAgICBzYW1wbGVzJGJldGEwICsgCiAgICBzYW1wbGVzJGJldGExICogZGZfY292aWQkaW5jb21lW2ldICsgCiAgICBzYW1wbGVzJGJldGEyICogZGZfY292aWQkYnNbaV0KICB5cmVwWywgaV0gPSBycG9pcyhuY3ljbGVzLCBsYW1iZGEgPSBleHAobG9nX2wpKQp9CgpwcGNfaW50ZXJ2YWxzKHkgPSB5LCB5cmVwID0geXJlcCwgeCA9IGRmX2NvdmlkJHN0YXRlX2FiYikgKyBnZ3RpdGxlKCJNb2RlbCBBOiBQb3N0ZXJpb3IgcHJlZGljdGl2ZSBpbnRlcnZhbHMgdG8gdGhlIG9ic2VydmVkIGRhdGEgdmFsdWVzIikgKyB4bGFiKCJTdGF0ZSIpICsgeWxhYigiQ292aWQgZGVhdGhzIikKCmBgYAoKCg==